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INTRODUCTION 



Shannon (1948) indicated how maximum entropy (ME) distributions can be derived by a 
straigtforward application of the calculus of variations technique. He defined the entropy 
of a probability density function p(x) as 



H = — p(x) \np(x) dx 



(1) 



Maximizing H subject to various side conditions is well-known in the literature as a 
method for deriving the forms of minimal information prior distributions; e.g. Jaynes 
(1968) and Zellner (1977). Jaynes (1982) has extensively analyzed examples in the 
discrete case, while in Lisman and Van Znylen (1972), Rao (1973) and Gokhale (1975), 
Kagan, Linjik continuous cases are considered. In the last case, the problem, in its 
general form, is the following 



maximize H — — I p(x) \np(x) dx 



subject to E{cj) n (x)} 



4> n (x)p(x) dx = fi n , n = 0,...,N 



(2) 



where /i — 1 , 4>o(x) = 1 an d 4> n (x),n — 0,...,N are N known functions, and fi n , n = 
0, . . . , N are the given expectation data. The classical solution of this problem is given 
by 



p(x) = exp 



N 



n=0 



(3) 



The (N + l) Lagrangien parameters A = [A , . . . , A„] are obtained by solving the follow- 
ing set of (N + l) nonlinear equations 



G n (X) 



[x) exp 



A' 



n=0 



dx = /!„, n — 0,...,N 



(4) 



The distributions defined by (3) form a great number of known distributions which are 
obtained by choosing the appropriate N and <f> n (x),n — 0, . . . , N. In general <f> n (x) are 
either the powers of x or the logarithm of x. See Mukhrejee and Hurst (1984), Zellner 
(1988), Mohammad-Djafari (1990) for many other examples and discussions. Special 
cases have been extensively analyzed and used by many authors. When <j> n (x) = x n ,n = 
0,...,N then /i„,n = 0,...,N are the given N moments of the distribution. See, for 
example, Zellner (1988) for a numerical implementation in the case N = 4. 

In this communication we propose three programs written in MATLAB to solve 
the system of equations (4). The first is a general program where <f> n (x) can be any 
functions. The second is a special case where <j> n {x) = x n ,n = 0,...,N. In this case 
the fi n are the geometrical moments of p(x). The third is a special case where (f) n (x) = 
exp(-jnux) , n — 0, . . . , N. In this case the fx n are the trigonometrical moments (Fourier 
components) of p(x). We give also some examples to illustrate the usefullness of these 
programs. 



PRINCIPLE OF THE METHOD 

We have seen that the solution of the standard ME problem is given by (3) in which the 
Lagrange multipliers A are obtained by solving the nonlinear equations (4). In general, 
these equations are solved by the standard Newton method which consists of expanding 
G n (X) in Taylor's series around trial values of the lambda's, drop the quadratic and 
higher order terms, and solve the resulting linear system iteratively. We give here the 
details of the numerical method that we implemented. When developing the G n (X) 
in equations (4) in first order Taylor's series around the trial A , the resulting linear 
equations are given by 

G n (X) - G n (A°) + (A - A )* [gradG' n (A)] (A=A o ) = pt n: n = 0,...,N (5) 
Noting the vectors S and v by 

<5 = A-A° 



OM* 



v=[ f i -Go(\°),...,fiN-G N (\ )} 



and the matrix G by 



G 



9nk 



dG n (X) 



n,k = 0,...,N 



(6) 



J (A=A°) 



then equations (5) become 



GS = v 



(7) 



This system is solved for 8 from which we drive A = A + 8, which becomes our new 
initial vector A and the iterations continue until 8 becomes appropriately small. Note 
that the matrix G is a symmetric one and we have 



9nk =9kn = - I <Pn{x) <Pk{x) exp 



N 



n=0 



dx n,k = 0,...,N (8) 



So in each iteration we have to calculate the N(N — l)/2 integrals in the equation (8). 
The algorithm of the general Maximum Entropy problem is then as follows: 

1 . Define the range and the discretization step of x ( xmin , xmax , dx ) . 

2. Write a function to calculate <p n (x) , n = 0, . . . , N ( f i n_x ) . 

3. Start the iterative procedure with an initial estimate A (lambdaO) . 

4. Calculate the (N + 1) integrals in equations (4) and the N(N — l)/2 distinct 
elements g nk of the matrix G by calculating the integrals in the equations(8) 

(Gn, gnk) . 

5. Solve the equation (7) to find 8 (delta) . 

6. Calculate A = A + <5 and go back to step 3 until 8 becomes negligible. 

The calculus of the integrals in equations (4) and (8) can be made by a univariate 
Simpson's method. We have used a very simplified version of this method. 



Case of geometrical moments 

Now consider the special case of moments problem where 4> n (x) 
0, . . . , N. In this case equations (3), (4) and (8) become 



x n , n 



p(x) = exp 



JV 



m=0 



X 



(9) 



G n (X) = J 



x exp 



N 



•5> 

m=0 



X 



dx = fi n , n — 0, . . . , N 



(10) 



9nk 9kn 





N 


j x n x k exp 






m=0 



dx 



-G n+k {\) n,k = 0,...,N (11) 



This means that the [(N + 1) x (JV + 1)] matrix G in equation (7) becomes a symmet- 
ric Hankel matrix which is entirely defined by 2N + 1 values G n (X),n = 0, . . . ,27V. So 
the algorithm in this case is the same as in the precedent one with two simplifications 

1. In step 2 we do not need to write a seperate function to calculate the functions 

<f) n (x)=x n ,n = 0,...,N. 



2. In step 4 the number of integral evaluations is reduced, because the elements g n k 
of the matrix G are related to the integrals G n (X) in equations (10). This matrix is 
defined entirely by only 2N + 1 components. 



Case of trigonometrical moments 

Another interesting special case is the case where the data are the Fourier components 

of p(x) 

E{exp (— jncuox)} = J exp (— jncoox) p(x) dx = fx n , n — 0,...,N, (12) 

where fi n may be complex-valued and has the property /i_ n = ii n . This means that we 
have the following relations 

(f) n (x) = exp(-jnu x), n = -N, . . . ,0, . . . N, (13) 



p(x) = exp 



N 



—Real A n exp (— jnuox) 



n=0 



(14) 



G n (X) = J exp (— jncu x) p(x) dx, n = 0, . . . ,N, (15) 

** - {-oZtw -*=° N - (16) 

so that all the elements of the matrix G are related to the discrete Fourier transforms of 
p(x). Note that G is a Hermitian Toeplitz matrix. 



EXAMPLES AND NUMERICAL EXPERIMENTS 

To illustrate the usefullness of the proposed programs we consider first the case of the 
Gamma distribution 

p(x;a,/3) = ^jz -x a exp(-[3x), x > Q,a < 1,(3 > 0. (17) 

r(l — a) 

This distribution can be considered as a ME distribution when the constraints are 

fp(x;a,P)dx =1 ( normalization 4>q(x) =1, 

f xp(x;a,[3)dx —Hi I 4>i(x) —x, (18) 

f\n(x)p(x;a,(3)dx — /i 2 { <f>2(x) =\n(x). 

This is easy to verify because the equation (12) can be written as 

p(x; a, (5) — exp [— A — Aix — A 2 ln(x)] 



0(l-a) 

with A = — In— r, Ai = (3 and \ 2 = —a. 

r(l — a) 

Now consider the following problem 

Given [i\ and // 2 determine A , Ai and A2. 

This can be done by the standard ME method. To do this, first we must define the range 
of x, (xmin, xmax, dx) , and write a function f in_x to calculate the functions 
4>o(x) = 1, 4>i(x) = x and fc^x) = hix (See the function f inl_x in Annex). Then we 
must define an initial estimate A for A and, finally, let the program works. 

The case of the Gamma distribution is interesting because there is an analytic relation 
between (a,f3) and the mean m = E{x} and variance a 2 = E{(x — m) 2 } which is 



or inversely 



m — (l — a)//3 
a 2 = (l-a)/(3 2 ' 



a = (a 2 — m 2 ) / a 2 
(3 = m/a 2 , 



(19) 



(20) 



so that we can use these relations to determine m and a 2 . Note also that the corre- 
sponding entropy of the final result is a byproduct of the function. Table (1) gives some 
numerical results obtained by ME_DENS1 program (See Annex). 



Table 1. 






a 


P 


m 


a 2 


0.2000 
0.2000 
0.3000 


-3.0000 
-2.0000 
-1.5000 


0.2156 
-0.4124 
-0.6969 


-3.0962 
-6.9968 
-5.3493 


0.2533 
0.2019 
0.3172 


0.0818 
0.0289 
0.0593 



The next example is the case of a quartic distribution 



p(x) = exp 



?1=0 



(21) 



This distribution can be considered as a ME distribution when the constraints are 

E{x n }= / x n p(x)dx — (j> n , ra = 0,..., 4 with fx — l. 



(22) 



Now consider the following problem : Given ii n ,n — 1, . . . ,4 calculate A n ,n = 0, . . . ,4 
. This can be done by the ME_DENS2 program. Table (2) gives some numerical results 
obtained by this program: 



Table 2. 



A*i yW2 At3 At4 A Ai 




0.2 0.05 0.10 0.1992 1.7599 
0.3 0.00 0.15 0.9392 0.000 
0.3 0.00 0.15 0.9392 0.000 



2.2229 -3.9375 0.4201 
-3.3414 0.0000 4.6875 
-3.3414 0.0000 4.6875 



These examples show how to use the proposed programs. A third example is also 
given in Annex which shows how to use the ME_DENS3 program which considers the 
case of trigonometric moments. 



In this paper we addressed first the class of ME distributions when the available data are a 
finite set of expectations \x n = E { 4> n (x) } of some known functions (f) n (x),n = 0,...,N. 
We proposed then three Matlab programs to solve this problem by a Newton-Raphson 
method in general case, in case of geometrical moments data where <f) n (x) = x n and 
in case of trigonometrical moments where (j) n (x) = exp(—jnu; x). Finally, we gave 
some numerical results for some special examples who show how to use the proposed 
programs. 
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ANNEX A 



unction [ lambda, p, entr ] =me_densl (mu, x, la: 
ME_DENS1 

[LAMBDA, P , ENTR] =ME_DENS1 (MU, X, LAMBDA ) 



This program calculates th 
probability density functi 
N contstraints in the form 
E{fin(x> }=MU{n) n=0:N 



is p(x) 



ipliers of the ME 
he knowledge of the 



Lth fiO(x)=l, MU(0)=1. 



straints MU(n) ,n=l:N. 

of the variation of x. 
st estimate of the LAMBDA s . 



MU is a table containing the 

g the r 
ing the 
optional. ) 

ing the resulting Lagrange parameters 
ing the resulting pdf p(x) . 
is a table containing the entropy values at each 
iteration . 



is a table contai 
is a table defini 
is a table contai 
(This argument 
is a table contai 
is a table contai 



Author : 
Date : 



A. Mohammad-D jafari 
10-01-1991 



21 mu=mu<:); 



; lx=length (x) ; 
xmin=x(l); xmax=x(lx); 



dx=x(2)-x(l); 



iflnargin == 2} 

lambda=zeros (size (mu) ) ; 

lambda ( 1 ) =log ( xmax-xmin 
else 

lambda=lambdaO ( : } ; 
end 

N=length (lambda) ; 



initialize LAMBDA 

% This produces 
distribution . 



33 fin=finl_x(x); 

34 % 

35 iter=0; 

36 while 1 

37 iter=iter+l; 

38 dispC 



f inl_x (x) is an external 
function which provides f 



start iterati 



' ) ; disp( [' iter=' , num2str (iter) ] ) ; 



p=exp (- (f i. 
plot <x,p> ; 



Calculate p (x) 
plot it 



G=zeros (N, 1) ; 
for n=l:N 

G(n)=dx*sum(fi 
end 



% Calculate Gn 



entr (iter) =lambda' *G (1 :N) ; 



% Calculate the entropy value 



disp ( [ ' Entropy=' , num2str (entr (iter) ) ] ) 



gnk=zeros (N, N) ; 

gnk (1, : ) =-G' ; gnk ( : , 1) =-G; 

for i=2:N 

for j=2:i 
gnk(i,j)=-dx*sum(fin(:,j). 

end 
end 

for i=2:N 



for j = 



+ 1 :N 



gnk(i, j)=gnk(j, 
end 
end 



% Calculate gnk 

% first line and 

% lower triangle 

% matrix G 



column 
of the 



triangle part of the 



=mu-G; 



Calculate v 
Calculate delta 
Calculate lambda 
Stopping rules 



delta=gnk\v; 
lambda=lambda+delta; 
eps=le-6; 

if (abs (delta . /lambda) <eps ) , break, end 

if (iter>2) 

if (abs ( (entr (iter) -entr (iter-1) ) /entr (iter) )<eps) , break, end 
end 



p=exp (- (fin* lambda) ) ; 

plot (x,p) ; 

entr=entr ( : ) ; 

dispC END 



Calculate the 
plot it 



1 % 

2 %ME1 

3 % This script shows how to use the function ME_DENS1 

4 % in the case of the Gamma distribution, (see Example 1.) 

5 xmin=0 . 0001; xmax=l; dx=0.01; % define the x axis 

6 x= [xmin : dx : xmax ] ' ; 

7 mu= [0 . 3, -1 . 5] ' ; % define the mu values 

8 [ lambda, p, entr] =me_densl (mu, x) ; 

9 alpha=-lambda(3) ; beta=lambda ( 2 } ; 

10 m= ( 1+alpha) /beta; sigma=m/beta; 

11 disp([mu' alpha beta m sigma entr (length (entr) )] } 

12 % 



1 function fin=finl_x (x) ; 

2 % This is the external function which calculates 

3 % the fin(x) in the special case of the Gamma distribution. 

4 % This is to be used with ME_densl. 

5 M=3; 

6 fin=zeros(length(x),M); 

7 fin(:,l)=ones(size(x)}; 

8 fin(:,2)=x; 

9 fin(:,3)=log(x); 
10 return 



function [lambda, p, entr ] =me_dens2 (mu, x, lambdaO ) 
%ME_DENS2 

% [LAMBDA, P , ENTR] =ME_DENS2 (MU, X, LAMBDAO) 

% This program calculates the Lagrange Multipliers of the ME 

% probability density functions p (x) from the knowledge of the 

% M moment contstraints in the form : 

% E{x A n)=mu(n> n=0:N with mu(0}=l. 

% 

i table containing the constraints MU (n ) , n=l : N . 
i table defining the range of the variation of x . 
i table containing the first estimate of the LAMBDAs . 
(This argument is optional . } 

a table containing the resulting Lagrange parameters, 
a table containing the resulting pdf p(x). 
a table containing the entropy values at each 
iteration . 



Author : 
Date : 



A. Mohammad-D jafari 
10-01-1991 



mu=mu ( : ) ; mu= [ 1 ; mu ] ; 
x=x(:}; lx=length(x) ; 
xmin=x (1) ; xmax=x (lx) ; 



% add mu(0}=l 
% x axis 
k(D; 



if(nargin == 2} 

lambda=zeros (size (mu) ) ; 

lambda ( 1 ) =log (xmax-xmin 
else 

lambda=lambdaO { : ) ; 
end 

N=length (lambda) ; 



initialize LAMBDA 

% This produces a uniform 
distribution . 



M=2*N-1; 

fin=zeros (length (x) ,M) ; 
fin( :, l)=ones (size(x) ) ; 
for n=2:M 



n)=x 



i(: 



-i); 



% Calcul de fin(x)=x 



iter=0; 

iter=iter+l; 
disp (' 



% start iterations 



' ) ; disp ( [' iter=' , num2str (iter) ] ) ; 



p=exp (- (f in ( : , 1 :N) 'lambda) ) ; 
plot (x,p) ; 



Calculate p (x) 
plot it 



G=zeros (M, 1) ; 
for n=l:M 

G(n)=dx*sum(fi 
end 



% Calculate Gn 



entr (iter ) =lambda' *G ( 1 :H) ; % Calculate the entropy value 
disp ( [ ' Entropy=' , num2str (entr (iter) ) ] ) 



gnk=zeros (N, N) ; 
for i=l:N 

gnk ( : , i ) =-G ( i : N+i-1 ) ; 
end 



% Calculate gnk 



G is a Hankel matr 



v=mu-G(l:N) ; 
delta=gnk\v; 
lambda=lambda+delta; 
eps=le-6; 
if (abs (delta . /lambda) <eps ) , 
if (iter>2) 

if (abs ( (entr (iter) -entr ( iter-1 ) ) /ent 
end 



Calculate v 
Calculate delta 
Calculate lambda 
Stopping rules 



break, end 
(iter) )<eps), break, end 



p=exp (- (fin ( : , 1 :N) 'lambda) ) ; 
plot (x,p) ; 
entr=entr ( : ) ; 

disp (' END ' ) 



% Calculate the 
% plot it 



%ME2 

% This sc 
% in the 
xmin=-l; x: 



x=l; dx=0.01; 



se the functi 
distribution 
% defin 



i ME_DENS2 
(see Example 2 . ) 
the x axis 



mu=[0.1, .3,0.1,-15]'; % 
[ lambda, p, entr] =me_dens2 (mu, x) ; 
disp ( [mu; lambda; entr (length (entr) ) ] ' ) 



define the mu values 



function [lambda, p, entr ] =me_dens3 (mu, x, lambdaO) 
%ME_DENS3 

% [ LAMBDA, P , ENTR ] =ME_DENS 3 (MU, X, LAMBDA ) 

% This program calculates the Lagrange Multipliers of the ME 
% probability density functions p (x) from the knowledge of the 
% Fourier moments values : 

% E{exp[-j n wO x])=mu(n) n=0:M with mu(0}=l. 



s a table containing the constraints MU(n),n=l:N. 
s a table defining the range of the variation of x. 
s a table containing the first estimate of the LAMBDAs . 
This argument is optional . ) 

s a table containing the resulting Lagrange parameters, 
s a table containing the resulting pdf p(x). 
s a table containing the entropy values at each 
iteration . 



Author 
Date 



A . Mohammad-D jaf ari 
10-01-1991 



mu=mu ( : ) ; mu= [ 1 ; mu ] ; 
x=x(:}; lx=length(x) ; 
xmin=x ( 1 ) ; xmax=x ( lx } ; dx 
if (nargin == 2} 

lambda=zeros (size (mu) ) ; 

lambda ( 1 ) =log (xmax-xmin) 
else 

lambda=lambdaO ( : } ; 
end 

N=length (lambda) ; 



% add mu (0) =1 
% x axis 
k(D; 

% initialize LAMBDA 

% This produces a uniform 
% distribution. 



32 M=2*N-1; 

33 fin=fin3_x(x,M); 

34 % 

35 iter=0; 

36 while 1 

37 iter=iter+l; 

38 dispC 



Calculate 
fin3_x(x) 



n(x)=exp[-jnw 
an external 



% function which provides fin(x) . 



% start iterations 



' ) ; disp ( [' iter=' , num2str (iter) ] ) ; 



p=exp (-real ( f i 
plot (x,p) ; 



% Calculate p (x) 
:N) ) 'real (lambda) +imag (f in ( : , 1:N) ) *imag (lambda) ) ; 
% plot it 



G=zeros (M, 1 ) ; % Calculate Gn 

for n=l:M 

G(n) =dx*sum(f in ( : , n) . *p) ; 
end 

%plot ( [real (G (1 : N) ) , real (mu) , imag (G (1 : N) ) , imag (mu) ] ) 

% 

entr (iter) =lambda' *G(1 :M) ; % Calculate the entropy 
disp ( [ ' Entropy=' , num2str (entr(iter) ) ] ) 



gnk= zeros (M, N) ; 
for n=l:M 

for k=l:n 
gnk (n, k) =-G (n-k+1) ; 

end 
end 



for 



=1 :N 



for k=n+l:N 

gnk (n, k) =-conj (G(k-n+l) ) ; 
end 
end 



Calculate gnk 

Matrix gnk is a Hermit i an 

Toeplitz matrix. 

Lower triangle part 



% Upper triangle part 



v=mu-G (1 : M) ; 
delta=gnk\v; 
lambda=lambda+delta; 
eps=le-3; 

if (abs (delta) . /abs (lambda) <eps) 



Calculate y 
Calculate delta 
Calculate lambda 
Stopping rules 



if (i 



>2> 



if (abs ( (entr (iter) -entr (iter 



break, end 
1) ) /entr (iter) ) <eps) , break, end 



end 

% Calculate p (x) 

p=exp(-real (fin( :, 1:N) ) *real (lambda) +imag (f in ( : , 
plot (x,p) ; % plot it 

entr=entr ( : ) ; 

disp (' END ' ) 

end 



1 :N) ) *imag (lambda) ) ; 



1 %ME3 

2 % This scripts shows how to use the function ME_DENS3 

3 % in the case of the trigonometric moments. 

4 clear; elf 

5 xmin=-5; xmax=5; dx=0.5; % define the x axis 

6 x= [xmin : dx : xmax] ' ; lx=length (x) ; 

7 p= (1/sqrt (2*pi) } *exp (- . 5* (x. *x) } ; % Gaussian distribution 

8 plot<x,p);title('p<x)'> 

9 % 

10 M=3;fin=fin3_x(x,M); % Calculate fin(x) 

11 % 

12 mu=zeros (M, 1 ) ; % Calculate mun 

13 for n=l:M 

14 mu (n) =dx*sum(fin ( :,n) . *p) ; 

15 end 

16 % 

17 w0=2*pi/ (xmax-xmin) ; w=w0* [0 :M-1] ' ; % Define the w axis 

18 % 

19 mu=mu<2:M); % Attention : mu(0) is added 

20 % in ME_DENS3 

21 [ lambda, p, entr] =me_dens3 (mu, x) ; 

22 disp ( [mu; lambda; entr ( length (entr) ) ] ' ) 



1 function f in=f in3_x (x, M) ; 

2 % This is the external function which calculates 

3 % the fin(x) in the special case of the Fourier moments. 

4 % This is to be used with ME_DEMS3 . 

5 % 

6 x=x<:}; lx=length (x) ; % x axis 

7 xmin=x(l); xmax=x(lx); dx=x < 2 } -x ( 1 ) ; 

8 % 

9 fin=zeros (lx,M) ; % 

10 fin{ :, l)=ones {size <x) } ; % fiO(x)=l 

11 w0=2*pi/ (xmax-xmin) ; jw0x= (sqrt (-1) *w0) *x; 

12 for n=2:M 

13 fin(:,n)=exp(-(n-l)*jw0x}; 

14 end 

15 return 



